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We perform full-potential screened-hybrid density-functional theory (DFT) calculations to com¬ 
pare the thermodynamic stability of neutral and charged states of the surface oxygen vacancy at 
the rutile TiO2(110) surface. Solid-state (QM/MM) embedded-cluster calculations are employed 
to account for the strong Ti02 polarization response to the charged defect states. Similar to the 
situation for the bulk O vacancy, the +2 charge state Vq'^ is found to be energetically by far most 
stable. Only for Fermi-level positions very close to the conduction band, small polarons may at 
best be trapped by the charged vacancy. The large decrease of the Vq'^ formation energy with 
decreasing Fermi-level position indicates strongly enhanced surface O vacancy concentrations for 
p-doped samples. 

PACS numbers: 71.15.Mb,68.47.Gh,68.55.Ln 


I. INTRODUCTION 

Its seemingly endless range of applications^'^ has made 
Ti02 one of the most studied transition metal oxides 
to date. Much of this material’s functionality in cor¬ 
responding (opto-)electronic, (photo-)catalytic or photo¬ 
voltaic applications derives not from the ideal bulk and 
surface structures, but is instead critically determined 
by intrinsic defects.^ Among these, the oxygen vacancy 
and in particular its nature as a charge trapping cen¬ 
ter have been most controversially discussed.^ The re¬ 
moval of an O atom from the bulk Ti02 lattice results 
in a single-particle defect state created from the three 
Ti dangling bonds that point into the vacancy. The en¬ 
ergetic position of the state depends sensitively on its 
electron occupancy and concomitant local lattice relax¬ 
ations. This occupancy can range from two electrons 
in the charge-neutral defect state (Vq), over one elec¬ 
tron in the singly-charged defect state (Vq), to empty 
in the doubly-charged state (Vq"^). Formerly prevalent 
seemingly contradictory schools of thought viewed the 
vacancy either as a shallow donor (Vq"^) that contributes 
to the n-type conductivity,or as an electrically inac¬ 
tive deep trap (Vq) that reduces neighboring Ti atoms 
and thus e.g. rationalizes an experimentally observed gap 
state about 0.8 eV below the conduction band.^^'^^ 

For bulk rutile Ti02 hybrid-functional density- 
functional theory (DFT) calculations by Janotti et 
^^ 15,16 r0(30]2tly resolved the preceding discrepancies by 
showing that Vq^ is the thermodynamically by far most 
stable configuration for all Fermi-level positions within 
the band gap. However, this charged donor state can 
trap one or two small polarons, in which excess electrons 
are localized on neighboring Ti^+ sites. The resulting 
weakly bound complex of shallow donor and small po¬ 
larons then naturally explains all experimental findings, 


but has to be carefully distinguished from the neutral 
Vq configuration in which the electrons occupy the de¬ 
fect state centered on the O vacancy site itself. 

As it is particularly surface O vacancies that play a cru¬ 
cial role in many of the Ti02 material applications^^'^^, 
it is important to assess how much of this novel un¬ 
derstanding derived for the bulk transfers also to these 
surface defects. Corresponding first-principles calcula¬ 
tions are, however, rather demanding. Already in the 
bulk case, at least hybrid-functional DFT is required to 
achieve an appropriate electron localization.Even 
with present-day computing power this constrains the 
system sizes that can be accessed. At the same time, 
the huge dielectric constant of Ti02 leads to a very large 
dielectric response in case of the charged defects. Here, 
not only lattice relaxations in the direct vicinity of the 
defect but the polarization of the entire semi-infinite sur¬ 
rounding medium contribute significantly. Within the 
conventional periodic boundary condition (PBC) super¬ 
cell approach this requires intricate extrapolation proce¬ 
dures involving supercells of increasing size.^^'^^ 

In this situation we instead opt for first-principles 
embedded-cluster calculations, in which the employed 
full-potential scheme allows for a numerically particu¬ 
larly efficient application of hybrid-functional DFT in¬ 
side the quantum mechanic (QM) cluster region.In 
the extended molecular mechanic (MM) embedding re¬ 
gion appropriately optimized interatomic potentials pro¬ 
vide a quantitative account of the strong Ti02 polar¬ 
ization response. We use this setup specifically to com¬ 
pute the formation energies, structural relaxations, and 
electronic structure of the bridging O vacancy at the ru¬ 
tile TiO2(110) surface. Consistent with the bulk calcula¬ 
tions of Janotti et al. we find that over a wide range of 
Fermi-level positions and oxygen chemical potentials the 
doubly-charged Vq"^ state is thermodynamically clearly 
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favored. The steeply decreasing formation energy of this 
charged defect with a lowering Fermi level then suggests 
p-doping as a promising avenue to tune the surface va¬ 
cancy concentrations and therewith the catalytic activity 
of this important material. 


II. METHODOLOGY 
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to 


too toQ 


4.587 

2.954 

111 257 

6.84 8.43 

HSE06^'> 

HSE06 (this work) 

4.588 

4.588 

2.951 

2.951 

278 402 
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MM (this work) 

4.587 

2.950 

3 

12 

5.76 6.73 


TABLE I. Rutile Ti02 lattice constants a and c, as well 
as its static (co, ^o) and high-frequency (e^ and e^) di¬ 
electric constants along the corresponding axes. Literature 
data from experiment^^“^^ and DFT-HSE06 calculations^^ are 
compared against our own calculations at the DET-HSE06 
level and with the parametrized interatomic potential for the 
MM region (see text). 
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EIG. 1. Schematic representation of the employed concen¬ 
tric solid-state embedding approach. A quantum-mechanical 
(QM) region hosting the surface oxygen vacancy is surrounded 
by a molecular mechanics (MM) region represented by a po¬ 
larizable interatomic potential. Spurious charge leakage out 
of the QM region is prevented through a transition shell in 
which cations are described with pseudopotentials (PPs). The 
full electrostatic potential of the infinite crystal surface is re¬ 
produced by placing point charges with fitted values around 
the MM region. 


We describe the localized surface defect in a concen¬ 
tric solid-state embedding approach as sketched in Fig. 1. 
The immediate vicinity of the defect and the rehybridiza¬ 
tion induced through it is treated at the QM level, and 
in particular through DFT. This QM region is embedded 
into a much larger MM region, which accounts for the 
longer-ranged dielectric properties of the TiO2(110) sur¬ 
face on the level of a polarizable interatomic potential. In 
a transition shell at the QM/MM boundary, cations are 
described by pseudopotentials (PPs) to prevent a spuri¬ 
ous overpolarization of the electron density (aka charge 
leakage). These PPs recover the long-range electrostatics 
of a point charge, but also have a repulsive short-range 
contribution which effectively mimics Pauli repulsion of 
core electrons. A final exterior shell of point charges is 
added at the outer boundary of the finite MM region. 
These point charges have values that are fitted to repro¬ 
duce the full electrostratic potential of the infinite surface 
inside the QM region.The next consecutive sections 
provide details of the DFT calculations, the parametriza- 
tion of the interatomic potential, and the employed sur¬ 
face models within this overall approach. 


A. Density-functional theory calculations 

All DFT calculations have been performed with 
the full-potential, all-electron framework 
Electronic exchange and correlation (xc) is treated at the 
level of screened hybrid-DFT, applying the HSE06 func¬ 
tional with the default mixing of 25% exact exchange 
and the default screening parameter of 0.2 Sys¬ 

tematic test calculations showed that the tier2 numerical 
atomic orbital basis set and the default tight settings for 
the atom-centered integration grids ensure a numerical 
convergence of the calculated defect formation energies 
within TlOmeV. Spin polarization is included through¬ 
out. As PPs for the transition shell we employ FHI98PP- 
generated^^ Ti^+ Kleinman-Bylander PPs^^ with non¬ 
local projector functions expanded up to the d-states. 
Further details of these potentials and their implementa¬ 
tion into FHI-aims can be found in ref. 25. 

The finite clusters describing the QM region are 
constructed using optimized bulk lattice parameters 
and relaxed positions of surface atoms at the ideal 
TiO2(110) surface as obtained from PBC-DFT super¬ 
cell calculations. The calculated lattice parameters of 
rutile Ti02 are listed in Table I and agree very well 
with experiments^ and preceding calculations at HSE06 
leveiss. The surface calculations employed symmetric 5 
0-Ti202-0 trilayer slabs, a c(4 x 2) surface unit-cell, and 
(4 X 4 X 1) Monkhorst-Pack k-point samplings^. Periodic 
slabs were separated in z-direction by 40 A of vacuum 
and electronically decoupled through a dipole correction. 
The surfaces were fully relaxed until residual forces fell 
below 10“S eV/A. 

Within PBCs the zero reference of the electrostatic 
potential is not uniquely defined.Notwithstanding, in 
surface calculations, the vacuum level can be determined 
as the electrostatic potential at the middle of the vac¬ 
uum separating the slabs. This gives access to the work 
function and the absolute position of the valence band 
maximum (VBM) at the surface, evBM(surf) = —8.2eV. 
To also access the bulk VBM position we computed the 
layer-resolved Tii^ core-level positions in an 11-trilayer 
thick (1x1) slab. These positions indicate a negligible 
band bending of the order of 30meV for the ideal stoi¬ 
chiometric TiO2(110) surface. This finding is confirmed 
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by additional calculations with up to 27-trilayer thick 
slabs using the PBE xc functional^^. For the purposes 
of this work we therefore equate surface and bulk VBM 
levels, and henceforth refer only to cvbm- 


B. Parametrization of the interatomic potential 


The HSE06 functional achieves a reliable account of 
the electronic contribution to the bulk dielectric prop¬ 
erties of rutile Ti02, as represented by the close match 
of the high-frequency dielectric constants and in 
Table L In contrast, it fails largely to describe the dom¬ 
inant lattice contribution additionally contained in the 
large static dielectric constants and e%. This has been 
traced back to its deficiencies in describing the intricate 
soft phonon modes of this material.In order to achieve 
a seamless embedding the employed interatomic poten¬ 
tials in the MM region should generally match the di¬ 
electric properties of the xc functional employed in the 
QM region. In the present case this would mean that 
the QM/MM approach then exhibits the same shortcom¬ 
ings with respect to the static dielectric properties; a 
point to which we will return in the discussion part be¬ 
low. In addition to the dielectric properties, the MM 
potential also has to match the QM lattice constants, to 
avoid artificially confining stress in particular when geo¬ 
metric relaxation of the QM region is to be considered. 
These demands highly challenge any existing interatomic 
potential.In the present case, this situation is further 
aggravated by the necessity to saturate the QM region 
with cationic norm-conserving PPs, which by definition 
have integer charges; in the case of Ti a charge of +4. 
For consistency, the remaining MM region is then also 
restricted to formal charges +4 and —2 on Ti cations 
and O anions, respectively. 

Oxygen ions in Ti02 are highly polarizable, and are in 
fact intrinsically polarized in the rutile structure. Using 
as interatomic potential a simple rigid ion model with 
formal charges does not capture this physics. In the 
QM/MM context, such potentials lead to an overpolar¬ 
ization at the QM cluster region boundary and an over¬ 
estimation of the electrostatic potential.^^’^^ In contrast, 
oxygen polarizability can be modeled efficiently within a 
polarizable shell-model^^, as has recently been demon¬ 
strated for Ti02 by Scanlon et Here, the oxy¬ 

gen anion is described by two point charges: A “core” 
charge representing the nuclei and closed-shell core elec¬ 
trons, and a “shell” charge simulating the valence elec¬ 
tron cloud. Mimicking electronic polarizability the oxy¬ 
gen core (c) and its shell (s) interact via a spring potential 
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where d^s is the distance between the core and shell 
charge, and kcs and r^s are parameters defining the 
potential. In this model^^, the dominant Coulomb inter¬ 
action between different oxygen shells (s-s) and between 


oxygen shells and Ti point charges (s-Ti) is furthermore 
augmented by Buckingham potentials 
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( 2 ) 


and 


A ( ds—Ti\ Cs—Ti 

Vs-Ti = A-Ti exp - - -g - , (3) 

V Ps-Ti J 4-Ti 

to provide an effective account of dispersive interactions 
and Pauli repulsion. Here, dg-s and dg-xi are the dis¬ 
tances between oxygen shells and between oxygen shell 
and Ti, respectively, and A, p and C are potential pa¬ 
rameters. Restricting the interaction between MM Ti 
cations to the mere formal charge electrostatics, the 
model is thus defined through a set of nine parameters: 
[^c—S5 '^c—S5 ^s—S5 Ps—si ^s—S5 ^s—Ti? Ps—Ti? Cg_xi: Qs\i with 
a final parameter for the charge on the oxygen shell. 

To obtain a seamless match to the QM region we per¬ 
form a global optimization of these parameters to rep¬ 
resent the bulk DFT-HSE06 lattice and dielectric con¬ 
stants. Specifically, we employ a differential evolutionary 
algorithm^^ from the Python package Inspyred 1.0^^ to 
minimize the dimensionless cost function 


rMM _ rDFT-HSE06 

_ 

^DFT-HSE06 

with Li = [a, c, e^, e^, e^, e^] and using the DFT-HSE06 
values from ref. 35, cf. Table I, as target values. In the 
corresponding MM calculations, the internal lattice pa¬ 
rameter u is always kept at its DFT-HSE06 value (0.305), 
and the static and high-frequency dielectric tensors are 
determined from the second derivative matrix of all MM 
particles or only of the shells, respectively.^^ 




kc-s [in eV/A 

Vc-s [in A] 

Qs [in e] 

23.67 

0.098 

-2.9332 

A [in eV] 

P [in A] 

C [in eV A“] 

s - s 23550 

0.2113 

38.55 

s- Ti 1838 

0.3207 

26.62 


TABLE II. Interatomic potential parameters, optimized to 
reproduce the bulk Ti02 DFT-HSE06 lattice parameters and 
high-frequency dielectric constants, see text. 

The best parameter sets generated this way still exhibit 
rather large errors, with F > 0.33. They typically ex¬ 
hibit substantial deviations in the lattice constants. This 
agrees not only with the observation from Cat low et alA 
who assigned the inability to reproduce both lattice and 
static dielectric constants to missing many-body terms 
in this class of interatomic potentials. It is also con¬ 
sistent with the finding of Lee et alA that deficiencies 
in the description of the Ti02 soft phonon modes (and 
therewith static dielectric constants) can be effectively 
cured through the use of different lattice constants. In 
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FIG. 2. Perspective view of the employed (a) Ti 22043 , (b) 
Ti32063 and (c) TueOgi clusters, each exhibiting a surface O 
vacancy in their central bridging O atom row. Ti atoms are 
shown as large white spheres, O atoms as small red spheres, 
and semi-transparent gray spheres mark the positions where 
PPs represent the immediately surrounding Ti-cations of the 
MM region. 


the present context, an accurate representation of the 
lattice constants is indispensable though, to avoid arti¬ 
ficial stress on the QM region. We therefore removed 
the static dielectric constants from the target set and 
immediately obtain parameter sets with significantly re¬ 
duced cost functions. The best parameter set exhibits 
an F = 0.007 and is used in all QM/MM calculations 
reported below. It is compiled in Table II and yields 
highly accurate lattice parameters and high-frequency di¬ 
electric constants as shown in Table I. Notwithstanding, 
its largely erroneous representation of the static dielec¬ 
tric properties (which is actually of the same magnitude 
but opposite sign to those obtained with HSE06 itself) 
is a concern and we discuss in Section III how this is 
addressed in our defect calculations. 


C. QM/MM setup 

In order to assess the convergence with respect to the 
employed QM region all calculations are done for a se¬ 
quence of three embedded cluster models, originally sug¬ 
gested by Ammal and Heyden^^. As shown in Fig. 2, all 
three clusters, Ti 22 043 , Ti 32063 and TUeOgi, are cen¬ 
tered on a bridging oxygen row, from where the central 
O atom has been removed to create the surface vacancy. 
Only the smallest and the largest cluster are of approx¬ 
imately hemispherical shape, while the middle Ti 32 063 
cluster is rather hemi-ellipsoidal, owing to the specifics 
of the rutile structure. Each cluster is embedded into 
an MM region that extends hemispherically up to a con¬ 
stant outer radius of 25 A. Eor the smallest cluster this 
translates to a total number of 3029 MM atoms, while 
for the largest cluster this translates to 2957 MM atoms. 
Every MM cation within 6 A vicinity of the QM clus¬ 
ter is replaced by PPs to suppress spurious charge leak¬ 
age out of the QM region, cf. Eig. 2. All MM cores are 
placed at positions according to those of the fully relaxed 
DET supercell reference calculation for the stoichiomet¬ 
ric TiO 2 ( 110 ) surface. All O shells are initialized to the 
fully relaxed state within a corresponding MM supercell 
reference calculation with the MM cores at exactly the 
same positions. A final exterior shell of 64 point charges 
around the MM region is then added, with charges fitted 
to reproduce the full electrostatic embedding potential 
of an infinite TiO2(110) surface.Einally, for every QM 
cluster all MM shells and QM atom positions are fully re¬ 
laxed. This setup defines what will henceforth be referred 
to as the ideal TiO2(110) surface. The surface defect se¬ 
tups are created from this reference setup, either without 
subsequent geometry relaxation (”non-relaxed” geome¬ 
tries) or with full geometry relaxation of all QM atoms 
except those at the QM region boundary until residual 
forces are below 10“^ eV/A (’’relaxed” geometries). 

Test calculations increasing the outer radius of the MM 
region up to 30 A show full convergence of the formation 
energy and electronic structure in the case of the neutral 
defect. In case of a net-charged QM region the polariza¬ 
tion response is much longer ranged though (oc q/{eR)). 
Eortunately, the missing polarization energy outside of 
the finite MM region can be reliably captured through 
an analytical correction. Eor a hemisphere with radius R 
in a continuum with a dielectric constant e and carrying 
a charge q at its center this analytic correction can be 
derived as^^ 

^ ( 5 ) 

Using the high-frequency dielectric constant in this ex¬ 
pression as further discussed below, we validated that 
with this post-correction also the formation energies of 
charged defects are fully converged with respect to the 
size of the MM region. Technically, we hereby use as 
isotropic dielectric constant the average over the diago¬ 
nal entries of the bulk dielectric tensor, e = (2e^+e^)/3. 
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As radius R we simply set the outer radius of the MM re¬ 
gion and variations of by ±1 A have a negligible effect 
on the calculated formation energies. As an approximate 
method we also use the analytical correction directly out¬ 
side the QM region {vide infra). Here, the choice of the 
radius is more critical and we determine it by measuring 
the semi-principal axes of the hemi-ellipsoid defined by 
the atomic positions of the QM cluster. R is then taken as 
the radius of a hemisphere with identical volume. Uncer¬ 
tainties of ±5 % in the thus determined radius translate 
in this case into an uncertainty of ±0.2 eV for A^^(Vq'^) 
and of ±0.05 eV for A^^(Vq). 

The actual QM/MM calculations for the thus de¬ 
fined setup are performed within the ChemShell 
environmentwith the interface described in detail 
before.We specifically use GULP^^ for the MM force 
calculations and the DL-FIND routine^^ for the geometry 
optimizations. Self-consistent polarization, aka shell op¬ 
timization, within the MM region as a response to an 
updated QM geometry is hereby calculated in a series of 
micro-iterations. 


D. Defect formation energies 

Neglecting vibrational entropic contributions in the 
solid state and defect-defect interactions in the dilute 
limit we define the formation energy to create a surface 
O vacancy in charge state q as^^ 

AEf(V^) = E(V^)-E(TiO2(110)) ± 

± Mo ± ± ^^poi(^) • (b) 

Here, E(TiO2(110)) and Biy^) are the total energies of 
ideal TiO 2 ( 110 ) and of TiO 2 ( 110 ) with the defect as ob¬ 
tained from our QM/MM calculations, respectively, mo 
is the chemical potential of oxygen and ep is the Fermi 
energy. In our sign convention a positive formation en¬ 
ergy implies a cost to create the defect. Correspondingly, 
the charge state exhibiting the lowest formation energy 
will be the thermodynamically stable state. 

Mo represents the energy of the reservoir which takes 
up the O atom that is removed from the crystal. It is gen¬ 
erally a variable. If the reservoir is a surrounding oxygen 
gas phase, mo is e.g. dependent on temperature and O 2 
pressure. Limits for /io can, however, be derived^^. In 
the extreme 0-rich limit 

Mo(O-rich) = 1/2F;(02) , (7) 

with E{ 02 ) the total energy of an isolated O 2 molecule. 
The opposite 0-poor (Ti-rich) limit can be assessed from 
the stability condition of bulk Ti 02 against decompo¬ 
sition into TUOs, Mo > —AF^^(Ti 203 ) ± 2 AF^^(Ti 02 ), 
where AF;f(Ti 203 ) and AF;f(Ti 02 ) are the formation 
energies of bulk Ti 203 and Ti02, respectively. At the 
HSE06 level this yields^^ 

Mo (O —poor) = /io(0 —rich) — 4.07 eV . (8) 


The Fermi energy ep is the energy of the reservoir 
where electrons released from the charged defects move 
to. Similar to mo? also ep is a variable that can e.g. be 
tuned through doping. It is most conveniently referenced 
with respect to evBM and below we thus report values 
from zero (Fermi level positioned at the VBM) to Aegap 
(Fermi level positioned at the conduction band minimum 
(CBM)), with Aegap the bulk band gap. In our QM/MM 
setup the zero reference for the Madelung potential and 
hence the VBM is by construction the vacuum level. In 
order to compute AE^ through eq. (6) the offset of the 
absolute ep value by evBM must therefore be considered. 

Within our supercell setup we calculate a Aegap = 
3.15 eV, which agrees well with preceding studies at the 
HSE06 level (3.19eV^^, 3.31 eV^^, 3.05eV^^, with the 
last study only employing 20% exact exchange). This 
also extends to the previously calculated absolute VBM 
position^^, which we determine at —8.2eV. The latter 
value is also consistent with the experimental value which 
can be estimated from the experimental work function 
and a presumed Eermi level position about 0.1-0.2eV be¬ 
low the CBM^^ in corresponding samples^^”^^. 

In an ideal world, evBM calculated for the different QM 
cluster sizes within our QM/MM setup would always be 
identical to the value obtained with the PBC supercell 
calculations. In practice, an imperfect QM/MM coupling 
affects the electronic structure and there in particular 
those orbitals with appreciable overlap with the QM/MM 
boundary. In the present case, this concerns precisely the 
highest occupied molecular orbital (HOMO), a.k.a. the 
VBM, of the QM clusters, which is of a rather delocal¬ 
ized nature. This results in a slight spurious electron 
depletion or aggregation in the QM region depending on 
the shape and size of the QM cluster, and therewith to 
shifts of the overall electrostatic potential. In order to 
compensate for this shift we apply a further correction 
term in the calculation of the absolute Eermi level posi¬ 
tion for eq. (6). This correction is determined once for 
each cluster size as the difference in the calculated un¬ 
relaxed, closed-shell singlet Vq defect-level position as 
compared to the corresponding level in the PBC calcu¬ 
lation. In contrast to the VBM/HOMO, the defect level 
is well localized in the center of the cluster. This min¬ 
imizes short-range artifacts from the QM/MM interface 
and probes (and corrects) the electrostatic potential ex¬ 
actly in the relevant region. 


III. RESULTS ^ DISCUSSION 
A. Electronic defect structure 

We begin our investigation by demonstrating the high- 
quality representation of the surface electronic struc¬ 
ture obtained with our QM/MM setup, which could 
generally not be achieved without the consistent re- 
parametrization of the MM potentials. Eigure 3 com¬ 
pares the density of states (DOS) for the ideal TiO 2 ( 110 ) 



6 



FIG. 3. Density of states (DOS) of the ideal TiO2(110) sur¬ 
face per Ti 02 formula unit as calculated with the embedded 
Ti 46 092 cluster (red) and with a supercell geometry (blue). 
evBM is used as zero reference and filled states are depicted 
in darker color. A Gaussian smearing (a = 0.1 eV) is applied. 


FIG. 4. Perspective view of the electron localization in the 
lowest-energy open-shell triplet conformation found for the 
neutral defect state with the embedded TueOgi cluster. Vi¬ 
sualized is the electron density contour at 0.01 e/A^. 


surface obtained with the embedded Ti46 092 cluster 
against the results from the slab reference calculation. 
Good agreement is achieved for both filled and empty 
states. This agreement extends not only to the band gap 
Acgap or to the width of the valence band. Both are 
reproduced within 0.5 eV and 0.9 eV, respectively. Also 
more subtle features within the bands are reflected rather 
well. Equivalent findings are obtained for the neutral 
defect, where a straightforward comparison to supercell 
calculations is also possible (not shown). 

For both the ideal surface and the neutral defect case 
the DOS calculated with the smaller QM clusters does 
also not differ much with respect to the corresponding 
calculation with the largest cluster. Intra-band features 
change, but the band gap or the width of the valence 
band varies each time by less than 0.2 eV. Even on abso¬ 
lute scale the spectra of the smallest and largest cluster 
match almost perfectly. In contrast, the intermediate 
Ti32 064 cluster reveals a potential offset by 0.4eV up 
in energy, which can be attributed to its prolate shape. 
This cluster exposes the largest surface to volume ratio 
to the QM/MM interface, and is hence most affected by 
imperfect QM/MM coupling. Note that this potential 
offset is compensated in eq. (6) through the correction 
procedure described in the previous chapter and, thus, 
does not affect the calculated formation energies. 

Introducing the oxygen vacancy gives rise to electronic 
defect states.Depending on the charge state of the de¬ 
fect these states are empty, singly or fully occupied. The 
two electrons in the fully occupied state of the Vq de¬ 
fect can thereby form three different electronic configura¬ 
tions with differing spin multiplicities: a closed-shell sin¬ 
glet, an open-shell singlet and an open-shell triplet. The 
closed-shell singlet corresponds to the paired excess elec¬ 
trons occupying the same orbital localized at the defect 
site. In the unpaired open-shell configurations both elec¬ 
trons occupy different orbitals of either parallel (triplet) 
or anti-parallel (singlet) spin. Rather than as intrinsic 


neutral Vq defect the latter configurations should thus be 
seen as a singly- (doubly-) charged defect with one (two) 
trapped small polaron(s).^^ In the limit of large electron 
separation both open-shell spin configurations must be¬ 
come energetically degenerate. Eor closer distances, the 
different electronic configurations give rise to different 
formation energies. In order to determine the appropri¬ 
ate neutral defect state energetics for comparison to the 
charged states, it is therefore necessary to identify the 
corresponding lowest-energy electronic configuration. 

As extensively discussed by Deskins et al the po- 
larons can be localized in numerous different patterns 
around the defect site with varying relative stabilities.^^ 
Applying the same initialization tricks as applied by De¬ 
skins et al to achieve these different localizations, we 
always end up with one electron trapped at the defect 
site, while the other electron can be trapped in different 
locations. Whenever these final sites are separated by 
more than 3 A the two different open-shell spin configu¬ 
rations emerge as energetically degenerate - exactly as 
reported by Deskins et al^^. In contrast to their obser¬ 
vation convergence to an open-shell singlet configuration 
proved essentially impossible for smaller electron separa¬ 
tions in our QM/MM setup. Even when starting from 
geometries optimized in the triplet state and with both 
excess electrons initially separated, the relaxed singlet 
geometries always converged to the closed-shell config¬ 
uration with paired electrons. Open-shell singlet states 
could only be stabilized in the largest QM cluster for 
rather large separations, then yielding energetically es¬ 
sentially degenerate stabilities to the open-shell triplet 
configuration as stated before. 

As further discussed below the closed-shell singlet con¬ 
figuration is energetically much less stable than any of 
the various open-shell localizations. In our case the lat¬ 
ter all exhibit rather similar relative stabilities within 
0.1-0.2eV. This is in disagreement to Deskins et 
and others^^’^^”^^, who reported relative stabilities vary¬ 
ing over a range of 1.76 eV and with the lowest-energy 
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conformations localizing both excess electrons in the sub¬ 
surface layer. Although in principle possible in the 
largest T^^Ogi QM cluster, we can not stabilize such 
conformations. In contrast, the lowest-energy conforma¬ 
tion found consistently across all three cluster sizes is 
depicted in Fig. 4, showing one polaron localized at the 
defect Ti-sites of the bridging row and one on a neighbor¬ 
ing five-fold coordinated surface Ti-atom with very little 
contributions at sub-surface Ti-sites. This difference to 
the preceding studies^^’^^”^^ could potentially arise out 
of remaining finite-size effects of our QM/MM setup at 
even the largest QM cluster employed. Alternatively, it 
could be due to their different treatment of the electronic 
structure at the effective GGA+U level. In this respect, 
we note that the surface localization site of our most 
stable conformation is in perfect agreement with previ¬ 
ous results by Di Valentin and co-workers employing the 
B3LYP hybrid functional. 

Having potentially missed the true lowest-energy solu¬ 
tion puts a certain uncertainty on the energetics of the 
open-shell neutral defect. However, as we will further dis¬ 
cuss below, even a higher stability by 0.87 eV that would 
then correspond to the sub-surface solution as reported 
by Deskins et al would not critically affect our conclu¬ 
sions on its relative stability with respect to the charged 
defects states. Even though therewith of no concern to 
the present study we nevertheless plan to further investi¬ 
gate this issue with large supercell hybrid functional cal¬ 
culations and GGA+U QM/MM calculations. We note 
that an important aspect are certainly the used lattice 
parameters, which in many GGA+U studies are simply 
those of the underlying GGA functional. We find the 
use of incorrect lattice parameters to critically affect the 
electronic level positions and therewith the polaron sta¬ 
bilities. Relying on the hybrid-functional optimized lat¬ 
tice parameters the electronic defect levels obtained in 
our QM/MM study are highly consistent with those of 
previous bulk and surface hybrid-functional studies. For 
all three QM clusters the closed-shell singlet state is lo¬ 
cated 0.8 eV below the GBM in the unrelaxed structure, 
and shifts up by 0.2 eV in energy upon including lattice 
relaxation, cf. refs. 15 and 21. The lowest-energy triplet 
solution exhibits two defect states 1.7 eV and 0.6 eV be¬ 
low the GBM in the unrelaxed structure. Upon lattice 
relaxation the electrons condense to two small polarons, 
the states of which at 1.1 and 0.9 eV below the GBM are 
fully consistent with experiment^’^^’^^'^^ {vide infra). 


B. Formation energies of unrelaxed vacancies 

In particular the large static dielectric response of Ti02 
is expected to significantly stabilize charged defects.In 
order to quantify this effect it is useful to first analyze 
the defect formation energies when structural relaxation 
in response to the defect is suppressed. From a method¬ 
ological point of view, this additionally provides the pos¬ 
sibility to assess an analytical scheme that may allow to 


Self-consistent MM electronic polarization 
Ti 22 043 Ti32063 TueOgi Slab 


Vq singlet 

VS 

v^+ 

5.60 

3.19 

3.98 

5.66 

3.22 

3.95 

5.60 5.59 

3.17 

3.95 


Analytic correction for electronic polarization 


Ti22 043 Ti32063 TUeOgi 

Vo singlet 

5.61 

5.66 

5.61 5.59 

vS 

3.19 

3.21 

3.17 

V?>+ 

3.98 

3.95 

3.95 


TABLE III. Unrelaxed surface defect formation energies (in 
eV) in the oxygen-rich limit and for a Fermi-level position 
at the VBM. Shown are results for the three different QM 
clusters employed, and, in the case of the neutral defect, also 
from the PBC supercell reference calculations. The data for 
the neutral Vq defect corresponds to the closed-shell singlet 
electronic configuration, which is the only one that can be sta¬ 
bilized in the absence of lattice relaxation. Upper rows corre¬ 
spond to a QM/MM setup, in which the electronic polariza¬ 
tion response to the defect is treated through self-consistent 
relaxation of the MM shells. Lower rows correspond to a 
QM/MM setup, where this response is approximately ob¬ 
tained through an analytical polarization correction approach 
(see text). 


compensate the error introduced by the inability of both 
the polarizable interatomic potential in QM/MM calcu¬ 
lations and the hybrid DFT functional in PBG supercell 
calculations to properly account for the large static Ti02 
dielectric constants, cf. Table 1. Gorrespondingly, Table 
HI summarizes calculated formation energies, in which 
thus only the high-frequency, electronic Ti02 polariza¬ 
tion is included that can accurately be accounted for by 
both the polarizable interatomic potential and the hy¬ 
brid functional. Shown is data for AE^ in the 0-rich 
limit and for the Fermi-level position at the VBM. The 
series of three QM cluster sizes demonstrates a fast con¬ 
vergence. For the neutral defect the obtained formation 
energy agrees furthermore perfectly with the correspond¬ 
ing value obtained by a PBG supercell calculation. 

In these QM/MM calculations the full electronic po¬ 
larization response of the material outside the QM re¬ 
gion is captured through the self-consistent relaxation of 
the O shells within the MM region, while the analytical 
correction AEpoi(g) in eq. (5) additionally accounts for 
the long-range contribution outside of the MM region. 
In these calculations the e value entering AEpoi(g) cor¬ 
responds - as appropriate for the unrelaxed case - to 
the isotropically averaged high-frequency dielectric con¬ 
stants, cf. Section HG. Table HI also compiles results ob¬ 
tained with a more approximate method, where the entire 
electronic polarization response of the material outside of 
the QM region is accounted for through the analytic cor¬ 
rection equation, eq. (5). Here, the O shells in the MM 
region are no longer allowed to relax after the creation of 
the defect. Instead, the analytical correction equation is 
employed with a radius R that corresponds to the outer 



radius of the QM region instead of the outer radius of 
the MM region used before. This method is numerically 
more advantageous as it does not need to achieve self- 
consistency between the MM shell polarization and the 
QM cluster. It also provides an elegant way to later on 
assess the error introduced by the erroneous description 
of the static dielectric constants through the polarizable 
force field. In principle, any value for e can be employed 
in A£’poi(g), i.e. also the ’’correct” experimental value. 
On the other hand, the method is more approximate and 
phenomenological. The polarization is only described at 
the isotropic continuum level and the outer boundary 
of the finite QM clusters is less well approximated by a 
hemisphere as the outer boundary of the much larger MM 
region. The results thus also depend more sensitively on 
the exact choice of the QM cluster radius R employed in 
A£’poi(g). For the data compiled in Table III the R for 
each cluster size has been chosen as a fit parameter to 
reproduce the corresponding self-consistent polarization 
results for both charged defects. The values obtained for 
R are for each cluster about 0.7-0.9 A larger than what 
would be obtained from the positions of the outermost 
QM atoms, cf. Section IIC, thus effectively accounting 
for the somewhat larger extension of the electron den¬ 
sity. Intriguingly, for each cluster size the same value 
of R can accurately reproduce the formation energies of 
both charged defects, cf. Table III. This shows that this 
approximate method can be applied without further sys¬ 
tematic errors and that for the polarization response out¬ 
side the QM region any electrostatic multipole moment 
of the electron density higher than the monopole term 
considered in AF^poi(g) can be indeed neglected. 

The data in Table III shows that already the compara¬ 
tively small electronic polarization is sufficient to largely 
stabilize the charged defects against the neutral one. We 
can quantify this stabilization through the calculated 
value of AEpo\{q) in case of the approximate method, 
i.e. when this term accounts for the entire response out¬ 
side of the QM region. Even for the largest QM cluster, 
where this additional stabilization through the far-range 
response is smallest, AEpoi(g) is still —0.60eV for the 
singly-charged Vq defect, whereas for the doubly-charged 
defect it even rises to —2.40eV. For the small¬ 
est QM cluster the corresponding values are —0.79eV 
and —3.16eV, respectively. For the neutral Vq defect 
AEpo\{q = 0) = 0. The essentially identical values for 
AE^(Vq) obtained with the self-consistent polarization 
and with the analytical correction method in Table III 
thus reveal that electronic stabilization outside the QM 
region is in this case negligible already for the smallest 
QM cluster. 

This different polarization response for neutral and 
charged defects is sufficient to make the latter thermo¬ 
dynamically more stable for a wide range of Fermi-level 
positions. The value of ep where the formation energies of 
charge state q and q' become equal define the transition 
level e^qjq'). With the values of Table III and eq. (6) we 
determine the transition level e(+/0) as 2.44 eV, i.e. less 


Self-consistent MM electronic polarization 
Ti22 043 Ti32063 TueOgi 


Vq singlet 

5.40 

5.50 

5.33 

Vq triplet 

4.89 

4.78 

4.83 

VS 

2.02 

1.83 

1.92 

v^+ 

0.43 

-0.20 

-0.37 


Analytic correction for full polarization 
Ti22 043 Ti32063 TUeOgi 


Vq singlet 

5.40 

5.50 

5.33 

Vq triplet 

4.89 

4.78 

4.83 

VS 

1.70 

1.56 

1.68 

V?>+ 

-0.69 

-1.27 

-1.30 


TABLE IV. Relaxed surface defect formation energies (in eV) 
for the three different QM clusters in the 0-rich limit and for 
a Fermi-level position at the VBM. Upper rows correspond 
to a QM/MM setup, in which only the electronic polarization 
response outside of the QM region is accounted for through 
self-consistent relaxation of the MM shells. Lower rows cor¬ 
respond to a QM/MM setup, where the full lattice and elec¬ 
tronic response is approximately obtained through an analyt¬ 
ical polarization correction approach (see text). 


than 1 eV below the CBM. Already when only account¬ 
ing for the (small) electronic response, the neutral defect 
level on which preceding calculations largely focused their 
attention is therefore not stable for Fermi-level positions 
over a large part of the band gap. When next also con¬ 
sidering the much larger lattice polarization of Ti 02 , this 
trend will be significantly enhanced, further disfavoring 
the neutral defect. 


C. Formation energies of relaxed vacancies 

The results of the preceding section already indicate 
the relevance of the charged defect states, and corre¬ 
spondingly the hitherto barely explored necessity to reli¬ 
ably describe them with electronic structure calculations. 
When moving to the fully relaxed defect formation ener¬ 
gies, our QM/MM scheme is limited by the shortcomings 
of the polarizable interatomic potential with respect to 
the static dielectric constants. Thus restricting struc¬ 
tural relaxation to the QM region (and maintaining the 
full electronic polarization as in the previous section), we 
obtain the defect formation energies compiled in Table 
IV. Comparing the entries for the neutral defect in Ta¬ 
ble III and IV we observe only a minute lowering of AE^ 
by 0.2 eV for the closed-shell singlet electronic configura¬ 
tion (henceforth denoted as singlet) at all three cluster 
sizes. The additional lattice response accounted for in 
the relaxed calculations is thus quickly converged over 
the few atomic shells contained in the finite QM clus¬ 
ters. While only a marginal quantitative effect for the 
singlet, it is only this short-range structural relaxation 
that stabilizes the neutral open-shell singlet and triplet 
configurations at all. Nevertheless, also for the degen¬ 
erate lowest-energy such configurations (henceforth de- 
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noted as triplet/singlet) quick convergence with the size 
of the QM region is achieved. 

In line with the expectations from the large static Ti02 
dielectric constants, the additional lattice-polarization 
stabilization of the charged defects is instead signifi¬ 
cant even when only accounting for the relaxation of the 
nearest-neighbor atomic shells in the QM region. Com¬ 
paring again the corresponding entries in Tables III and 
IV, this stabilization is of the order of 2 eV for the singly- 
charged defect and amounts even to around 4 eV for the 
doubly-charged defect. As polarization of the environ¬ 
ment scales with the square of the charge, cf. eq. (5), 
this drastic increase is not unexpected. At the latest 
for the Vq^ defect no convergence can concomitantly be 
achieved over the present range of QM cluster sizes. In 
view of the strong lattice relaxations calculated for this 
defect, this is also not surprising. With an outwards re¬ 
laxation of up to 0.45 A the nearest-neighbor Ti atoms 
around the defect exhibit the largest displacements. Cor¬ 
responding displacements in the two larger Ti32063 and 
TUeOgi QM clusters differ by less than 0.05 A, proving 
that this shortest-range polarization contribution is well 
converged at these cluster sizes. However, in the largest 
Ti4609i QM cluster the next two shells of Ti atoms still 
show maximum displacements of 0.12 A and 0.05 A, re¬ 
spectively. These non-negligible relaxation contributions 
can no longer be captured with the smaller clusters. 

Just accounting for the lattice relaxation inside the fi¬ 
nite QM region is therefore not sufficient to reliably deter¬ 
mine the formation energies of the charged defects. The 
relaxation patterns possible inside the largest TUeOgi 
QM cluster would require a c(2 x 5) surface unit-cell in 
PBC supercell calculations. Due to the periodic images 
the relaxation patterns outside this unit cell are spu¬ 
rious. Our findings for the QM/MM setup are there¬ 
fore paralleled by the equivalent insight that even cor¬ 
respondingly large surface unit-cells are not sufficient to 
absolutely converge charged defects in PBC calculations. 
As these system sizes are at the upper limit of what is 
presently tractable at the hybrid functional level, extrap¬ 
olation procedures are thus required within the PBC ap¬ 
proach. Within our QM/MM approach we can find an 
analogue in the analytic correction procedure described 
in the last section. Indeed, using exactly the same ra¬ 
dius R for each of the three cluster sizes as established 
in the unrelaxed calculations we can reproduce the 
of Table IV with the same accuracy as was the case in 
Table III (not shown). Here, this means that we used for 
e in AEpo\{q) of eq. (5) exactly the high-frequency value 
that corresponds on average to the ones of the polarizable 
interatomic potential. 

With this confidence we apply the approximate an¬ 
alytic correction also for the calculation of fully relaxed 
formation energies. For this, we now use as e in AEpo\{q) 
the isotropic average of the large static dielectric con¬ 
stants derived from HSE06, cf. Table I. Keeping exactly 
the same radii R as before in the calculation of AEpoi(g), 
this approach thus accounts for the full lattice and elec¬ 


tronic relaxation explicitly in the QM region, and ac¬ 
counts for the same full lattice and electronic polarization 
on the continuum level outside of the QM region. The re¬ 
sults compiled in Table IV indicate another sizable stabi¬ 
lization in particular for the doubly-charged Vq^ defect, 
the formation energy of which is now also converged over 
the two larger QM clusters. Likely, the obtained AE^ 
values nevertheless still slightly overestimate the true re¬ 
laxed formation energies, as the outermost atoms in the 
QM region are not allowed to relax to avoid artifacts from 
an imperfect QM/MM coupling. From small variations 
of the employed radius R in AEpo\{q) to also effectively 
account for the corresponding shell of atoms, we estimate 
that this relaxation would lower the formation energy of 
Vq (Kq^) by another 0.1 eV (0.3eV). 

We therefore arrive at final values for the formation 
energies in the 0-rich limit and for a Fermi-level position 
at the VBM of 5.3 ± 0.1 eV (Vq singlet), 4.8 ± 0.1 eV 
(Vq triplet/singlet), 1.6±0.2eV (Vq) and —1.6±0.3eV 
(Vq"^). The stated error bars hereby reflect conserva¬ 
tive estimates accounting for the uncertainties implied 
by the analytic correction approach and the convergence 
with the QM cluster size. This uncertainty in particular 
for the doubly-charged vacancy is not entirely satisfying. 
Notwithstanding, we note that a similar uncertainty did 
arise in the calculation of the bulk defect formation en¬ 
ergies through the extrapolation procedure required in 
the employed PBC supercell approach^^, whereas a re¬ 
liable polarization-converged calculation of AE^ of the 
charged surface defects at the hybrid-functional level has 
never even been attempted. Intriguingly, our results are 
on the contrary very robust against the inaccuracy in¬ 
troduced by the still large deviation of the static HSE06 
dielectric constants with respect to experiment, cf. Ta¬ 
ble I. As apparent from the nature of eq. (5) values 
for e smaller than the HSE06-derived value will cause 
large variations in the long-range polarization correction. 
Varying the e employed in the analytic correction from 
the small value corresponding to the MM high-frequency 
dielectric constants to the one corresponding to the static 
HSE06 dielectric constants thus led to the just discussed 
large changes of the formation energies of the charged de¬ 
fects. However, a further increase of e to match the exper¬ 
imental dielectric properties of Table I instead yields only 
a negligible further stabilization of the singly- (doubly-) 
charged defect by 0.01 eV (0.03eV). 

In total, structural relaxation therefore lowers the for¬ 
mation energies of the closed-shell Vq singlet by 0.3 eV, 
of Vq by 1.6 eV, and of Vq"^ by 5.6 eV. In the work by 
Janotti et al for the bulk O vacancy the correspond¬ 
ing values were 0.3 eV, 0.9 eV and 3.5 eV.^^ This sug¬ 
gests a significantly stronger stabilization of the charged 
defects at the surface, which one may attribute to a gen¬ 
erally larger structural flexibility of surface atoms. In 
light of the preceding discussion on the insensitivity of 
the polarization-response to small variations away from 
the HSE06 dielectric properties, this conclusion should 
not be affected by the use of a tailored HSE functional 
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with a slightly different exact exchange mixing in this 
preceding work. More likely, the different functional will 
affect the description of the short-range QM rehybridiza¬ 
tion and therewith prevent a direct comparison of the 
absolute defect formation energies in the bulk and at the 
surface. Compared to our surface values of 5.3 eV (Vq 
singlet), 1.6eV (Vq) and —1.6,eV (Vq"^), Janotti et al. 
report 5.3eV, 2.2eV and —1.4eV for the bulk, respec¬ 
tively.^^ Overall, the numbers are intriguingly similar. 
Considering also the uncertainties due to PBC extrapola¬ 
tion or QM/MM convergence we can, however, not make 
any more detailed statements. Even though the calcu¬ 
lated AE^(Vq'^) reffect the same trend we can therefore 
not comment on the strong surface segregation tendency 
obtained previously at the GGA or GGA+U level. 
Corresponding (polarization-response unconverged) cal¬ 
culations had predicted a higher stability of the surface 
O vacancy by ^ 0.7 — 1 eV as compared to its counterpart 
in the bulk. 


D. Thermodynamic stability 
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FIG. 5. Calculated relaxed surface defect formation energies 
as a function of the Fermi-level position in the band gap. 
evBM is used as zero reference for ep- The y-axis scales on 
the left and right correspond to AE^ in the the 0-poor and 
0 -rich limit, respectively, and are correspondingly shifted by 
4.07 eV, cf. eq. (8). 


As clear from eq. (6) the defect formation energies de¬ 
pend on the reservoirs available for the removed species, 
i.e. the oxygen chemical potential and the electronic 
Fermi level. Figure 5 therefore extends the previous dis¬ 
cussion for the 0-rich limit and for a Fermi-level position 
at the VBM, and shows the different AE^ as a function 
of the Fermi-level position and with ^-axis scales for both 


the 0-rich and 0-poor limit. For any Fermi-level position 
within the band gap, the doubly-charged defect ex¬ 

hibits the lowest formation energies; a situation that was 
equivalently obtained before for the bulk O vacancy. 
When also accounting for the full lattice relaxation, the 
transition levels e(2 + /+) and e(2 + /O) are thus located 
above the CBM, i.e. the surface O vacancy is a shallow 
donor. 

Even in the limit of strong n-type doping with a 
concomitant Eermi level close to the CBM the ener¬ 
getic gap to the neutral closed-shell singlet Vq with two 
electrons bound in the localized defect orbitals is still 
quite large. In full agreement with previous theoret¬ 
ical studies^^’^^’^^”^^ the equally charge-neutral open- 
shell triplet/singlet is obtained as an energetically pre¬ 
ferred electronic configuration. However, in contrast to 
the intrinsically neutral defect, this would instead cor¬ 
respond to a situation where the charged defect has 
trapped two small polarons^^’^^. In our calculations, this 
situation only becomes energetically degenerate to the 
bare doubly-charged Vq"^ defect for a Eermi-level posi¬ 
tion right at the CBM, cf. Eig. 5. Here, we have to recall 
though that Deskins et al reported a polaron localiza¬ 
tion in the sub-surface layer, which we could not con¬ 
firm with the present hybrid-functional QM/MM setup 
{vide supra). In their scheme, this sub-surface configura¬ 
tion was 0.87 eV more stable than the surface polaron 
triplet/singlet configuration we obtain as most stable, 
cf. Eig. 4. If the sub-surface solution by Deskins et al. 
is indeed physical, this would thus lower the open-shell 
triplet/singlet line in Eig. 5. The resulting lowering of 
the transition level e(2 -1- /O) to 0.44 eV below the CBM 
would then indicate the possibility of polaron trapping 
at the surface defects for corresponding Eermi-level posi¬ 
tions. As discussed by Janotti et al.^^ it is such trapped 
polarons (trapped possibly at the surface O defect but 
equally at other surface and bulk defects), not the O 
surface defect itself that are responsible for the defect 
state in the band gap observed in numerous experimen¬ 
tal studies^^. 

Eor Eermi-level positions further away from the CBM, 
the formation energy of the Vq"^ defect decreases rapidly, 
cf. Eig. 5. Under 0-rich conditions it is only this steep 
lowering of AE^ that eventually leads to values consistent 
with appreciable defect concentrations. This is consistent 
with experimental reports on increased defect concentra¬ 
tions upon p-doping.^^“^^ Equivalent reductions of AE^ 
(and concomitant increases of the defect concentration) 
can, of course, also be achieved with less 0-rich condi¬ 
tions, i.e. by lowering the O chemical potential. A strong 
presence of bridging-oxygen vacancies at the TiO 2 ( 110 ) 
surface has indeed frequently been observed in ultra-high 

vacuum experiments.apparent from Eig. 5 the 
lowering of AE^ with lower ep or /io is in fact so strong, 
that we eventually obtain negative formation energies. 
This unphysical result indicating a lattice instability is 
an artifact of the persistent use of defect formation en¬ 
ergies calculated for the dilute limit. Under conditions 
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corresponding to increased defect concentrations, defect- 
defect interactions as well as the build-up of a space- 
charge region^^ would in reality modify the formation 
energies to suppress negative values/^ In view of the 
huge energetic preference of the doubly-charged defect 
presently obtained for such conditions, it is, however, 
unlikely that this will change the energetic ordering of 
neutral and charged defects. 

IV. SUMMARY AND CONCLUSIONS 

We presented a solid-state QM/MM approach with 
a polarizable interatomic potential optimized to match 
the DFT xc functional employed in the QM region and 
analytical corrections for long-range polarization effects. 
Corresponding embedded-cluster calculations provide a 
first determination of the defect formation energies of 
neutral and charged O vacancies at the TiO2(110) sur¬ 
face at the hybrid-functional DFT level and containing a 
converged contribution of the strong dielectric response 
of this material. The stabilization of the singly- (V^(^) 
and doubly- (V^q^) charged surface defect through lat¬ 
tice relaxation is indeed found to be sizable, i.e. of the 
order of 1.6 eV and 5.6 eV, respectively. It is thus even 
larger than previously obtained for the bulk O vacancy^^, 
a fact that we attribute to the generally larger structural 
flexibility of surface atoms. 

The stabilization of in particular the doubly-charged 
defect is large enough to make it the thermodynam¬ 
ically most stable state for any Fermi-level position in 
the band gap. However, under the uncertainties of our 
approach we cannot exclude a possible trapping of small 
polarons at the charged defect for Fermi-level positions 


close to the CBM. The situation for the surface O va¬ 
cancy would then be fully equivalent to the one discussed 
by Janotti et al before for the bulk.^^’^^ The surface O 
vacancy is thus a shallow donor and the electronic de¬ 
fect state within the band gap observed experimentally 
results from trapped polarons, not from the intrinsic O 
defect itself. 

Within the nature of a charged defect, the formation 
energy of the surface O vacancy varies with the 

Fermi-level position in the band gap. In line with 
experimental reports this predicts largely increased 
vacancy concentrations upon p-doping. As surface 
vacancies are frequently discussed as reactive centers, 
systematic variations of doping concentrations may 
therefore provide an important avenue to tune the cat¬ 
alytic activity of Ti02. The presented embedded-cluster 
approach allows to efficiently address charged reaction 
intermediates or their binding to charged defects both 
under a converged account of the large polarization 
response and at hybrid-functional DFT or even beyond. 
This predestines the approach to quantitatively assess 
this avenue for this and other (polarizable) metal oxides. 
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